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Motivated by the question of how to pattern a surface in order to best speed nucleation from 
solution, we build on the work of Page and Sear [Phys. Rev. Lett. 97, 65701 (2006)] and calculate 
rates and free energy profiles for nucleation in the 3d Ising model in the presence of cuboidal pores. 
Pores of well-chosen aspect ratio can dramatically speed nucleation relative to a planar surface made 
of the same material, while badly-chosen pores provide no such enhancement. For a given pore, the 
maximum nucleation rate is achieved when one of its two horizontal dimensions attains a critical 
length, largely irrespective of the other dimension (provided that the latter is large enough). This 
observation implies that patterning a surface in a raster-like fashion is a better strategy for speeding 
nucleation than e.g. scoring long grooves in it. 



I. INTRODUCTION 

The presence of a pore or a pit in a surface can affect 
nucleation in a profound way [TH5]. Page and Sear [9] 
used the 2d Ising model to demonstrate that nucleation 
of a new phase can be much faster in a pore than on a 
fiat surface of the same material, because nucleation can 
start at energetically-preferred binding sites in a pore 
corner. Further, for given thermodynamic conditions, 
they showed that there exists a pore width that maxi- 
mizes nucleation rate. The existence of this maximum 
follows immediately from the fact that as one makes a 
2d pore wider, the rate for nucleation into the pore is 
reduced, while the rate for nucleation out of the pore 
(into solution) is enhanced. Here we show that similar 
arguments in three dimensions suggest simple strategies 
for patterning a surface in order to best speed nucleation 
from solution. 

In making this claim, we present an analysis of nucle- 
ation in the Ising model that complements several previ- 
ous studies. We calculate rates and free energy profiles 
for nucleation in the 2d Ising model in the bulk, at fiat 
surfaces, and in the presence of rectangular pores, and 
in the 3d Ising model in the presence of cuboidal pores. 
We first calibrate our free energy sampling procedure 
by following Ref. [10] and comparing free energy pro- 
files for bulk nucleation in the 2d Ising model with the 
predictions of classical nucleation theory (CNT) mod- 
ified to accommodate nucleus shape fluctuations. As 
did Ref. [10 , we find excellent agreement between the- 
ory and simulation over a wide range of conditions. We 
next show that a flat surface can enhance nucleation, 
provided that this surface exerts a sufficiently large at- 
traction (explicit or effective) for the nucleating phase. 
Otherwise, nucleation happens in the bulk. While intu- 
itively reasonable, and described qualitatively by CNT, 
we argue that the need for such an attraction is obscured 
in the spin-spin representation of the Ising model, and 
is made clear only in the lattice gas (particle-vacancy) 
one. We then revisit the study of Ref. [9 by calculat- 
ing free energy profiles for pore nucleation in the 2d 
Ising model. These profiles support and complement 
the nucleation rates presented in that work, showing 
the existence of a pore size optimum for speeding nu- 



cleation from solution. We end with our main results, 
rates and free energy landscapes for nucleation in the 
3d Ising model in the presence of a cuboidal pore. We 
find that nucleation profiles display single or double bar- 
riers, depending on pore size and aspect ratio. Pores 
of well-chosen aspect ratio can dramatically speed nu- 
cleation relative to a planar surface made of the same 
material, while badly-chosen pores provide no such en- 
hancement. For a given pore, the minimum barrier to 
nucleation is achieved when one of its two horizontal 
dimensions attains a critical length, largely irrespective 
of the other dimension, provided that the latter is large 
enough. This observation implies that patterning a sur- 
face in a raster-like fashion is a better strategy for speed- 
ing nucleation than e.g. scoring long grooves in it. 



We note at the outset that there are important lim- 
itations to our study in terms of its relevance to real 
systems. We have chosen to work with the Ising model, 
a prototypical description of phase change and nucle- 
ation [TTHT6] , because it captures important effects of 
fluctuations and geometry crucial to many physical pro- 
cesses. Insights derived from it are often transferrable 
to many different physical systems [T7J [18] . In partic- 
ular, the result that there exists a pore size and shape 
in three dimensions that minimizes free energy barri- 
ers to nucleation follows immediately from geometrical 
considerations, in much the same way as the result that 
in two- and three-dimensional bulk space there exists a 
free energy barrier to nucleation. We expect therefore 
that this result should be relevant to three-dimensional 
systems generally. However, we do not represent impor- 
tant physical processes that may, in real systems, act 
to mask the existence of such an optimum. By repre- 
senting the nucleating phase as a lattice-based, struc- 
tureless one, we cannot capture potentially important 
effects like the mismatch in registry between a crystal 
and its template [2]. Furthermore, our simulation pro- 
tocol (grand-canonical Monte Carlo) is an efficient way 
of mapping free energy profiles, but ignores effects of 
mass transport and particle correlations that may be 
crucially important near real pores and surfaces. With 
these caveats in mind, we proceed to our study. 
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II. MODEL AND SIMULATION METHODS 

We consider homogeneous and heterogeneous nucle- 
ation in the 2- and 3-dimensional Ising model on a 
square or cubic lattice. Because we have in mind the 
nucleation of particles from solution, we find it conve- 
nient to work in the lattice gas (particle- vacancy) repre- 
sentation. Regardless of dimension, the energy function 
of our system is 

wall 

E = -J^UiUj - fjb^^m- Js^riinJ. (1) 

(ij) i ij 

Here rii = if site i is vacant, and rii = 1 if site i is 
occupied by a particle. The first two terms are the usual 
bulk ones: J is the strength of the nearest-neighbor 
coupling, and \i is a chemical potential that can be tuned 
to favor particles or vacancies. The first sum runs over 
all distinct nearest-neighbor bulk bonds, and the second 
sum runs over all bulk sites. The third term describes 
interactions between particles and walls: this sum runs 
over all bonds connecting wall sites to bulk sites. Wall 
sites are considered to be particles, i.e. nj = 1. Our 
simulations were done in the lattice gas representation, 
but we will use regular Ising (spin-spin) variables where 
convenient. For reference: via the usual mapping, rii = 
7,(1 + Si), where Si = ±1, the bulk lattice gas maps 
(ignoring constant terms) to the bulk Ising model 

E^^-Kj^SiSj-hJ^Si, (2) 

(ij) i 

where J = AK and /i = 2h — 2zK. Here z = 2d is 
the coordination number of the d-dimensional square or 
cubic lattice. 

We carried out simulations using a standard grand 
canonical Metropolis Monte Carlo (MC) procedure [19]. 
Each trial move consisted of an attempted change of 
state of a randomly-chosen lattice site. Trial moves re- 
sulting in an energy change AE were accepted with 
probability min (1, exp(— /3AE)), where f3 = l/(k&T). 
To study homogeneous nucleation in 2d, we used a sys- 
tem of iVbox = 100 2 sites with periodic boundary con- 
ditions in both directions. To model a wall we use pe- 
riodic boundaries in the horizontal direction only, and 
created a closed boundary along the top and bottom of 
the box by filling all sites in the first row of the lat- 
tice with particles. Following the protocol of Ref. [9], 
we also simulated rectangular pores in a system of 60 2 
sites. The depth of a pore was fixed at 30 sites and 
its width was allowed to vary. An analogous protocol 
was used in 3d simulations to generate cuboidal pores 
of fixed depth and a range of different length-to-width 
aspect ratios. 

To calculate free energy landscapes for nucleation we 
used standard umbrella sampling [20] protocols, and for 
the calculation of nucleation rates we used forward flux 



sampling (FFS) [21]. Brief details of both methods are 
given in Appendix | A 1| 



III. RESULTS: 2D SIMULATIONS 

A. Free energy barriers for bulk nucleation in the 
2d Ising model are well described by (modified) 
classical nucleation theory 

By way of calibration, it is instructive to compare 
nucleation free energy barriers calculated by umbrella 
sampling with the predictions of classical nucleation the- 
ory [22] [23] (CNT). Ref. [10, showed that sampled free 
energy barriers in the 2d Ising model are in excellent 
agreement with the predictions of the CNT-like expres- 
sion 

G C nt(N) = -AgN + 2aVTrN + Tk B TlnN + d(T). (3) 

Here G is the excess free energy of a droplet of N 'up' 
spins in a background of 'down' spins. The first term 
of this expression is the conventional bulk reward for 
growth of a circular droplet. Ag is the bulk free en- 
ergy difference between the two bulk phases, equal to 
2h at low temperature. At higher temperature (still 
below the critical one) the viable bulk phases are less 
dense than the all-up and all-down spin limits. Here 
we expect Ag ^ hAm to be a reasonable approxima- 
tion, where Am is the magnetization difference between 
bulk phases. For the conditions considered in this sec- 
tion, these estimates differed by at most about 3%; we 
therefore set Ag = 2h. The second term of Eq. ([3| 
is the surface tension penalty for growth of a circular 
droplet, a is the inter-phase surface tension. In 2d, 
the Ising model surface tension (at h = 0) in the direc- 
tion of either lattice vector is known from the Onsager 
solution [24], and is 

cry =2K- k B T In coth(f3K). (4) 

Because a non-square droplet cannot be accommodated 
perfectly on a square lattice, it is also useful to con- 
sider the orientationally-averaged effective droplet sur- 
face tension a e ^(T) derived by Shneidman et al. [25] : 

a eS (T)^^=(cr\\+a diag ), T>0.25T C . (5) 
Here 

X (T)= (l-sinh- 4 (2/3X)) 1/8 , (6) 

and 

c^diag = V2k B T lnsinh(2/3iT) (7) 

is the surface tension in the direction of the unit cell di- 
agonal [26] . In what follows we compare our simulations 
with the predictions of Eq. ^ using a defined both by 
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FIG. 1: Free energy barriers for homogeneous nucleation in the 2d Ising model are well-described by Eq. (pi). We start from 
an empty box (a box of vacancies), and use umbrella sampling to grow clusters of connected particles (droplets). Panels 
(a) and (c) show droplet free energy profiles as function of droplet size N, for various values of /cbT and for magnetic field 
strengths ft = 0.1 and ft = 0.06. Open circles denote the results of umbrella sampling simulations; solid lines are obtained 
from Eq. using the Shneidman et al. effective surface tension. In agreement with Ref. [10], the theoretical predictions 
match the simulation data closely. Panels (b) and (d) compare free energy barriers obtained from simulation in the plots 
directly above them. The dotted lines show the conventional (uncorrected) CNT barrier height predictions of Eq. ^ with 
t — d — 0. In all cases the uncorrected theory significantly underestimates the barrier height. Upward pointing triangles 
denote the results of Eq. ^ using the Onsager solution to the surface tension. Agreement is good at high temperature but 
less good at low temperature, where droplets are anisotropic (see snapshots on diagram). 



Eqns. Q and 

The third term of Eq. ([3| accounts for shape fluc- 
tuations of the droplet, and can be derived from field 
theoretic considerations of nucleation rates [27U30] . The 
shape fluctuation parameter r = 5/4 in 2d [51] . One im- 
portant contribution of Ref. [TO] was to recognize that 
this term can be considered a contribution to the free 
energy of a droplet. Without such a contribution, CNT 
and umbrella sampling are in quantitative disagreement. 
The final term of Eq. ([3| accounts for the fact that the 
conventional CNT expression has no clear origin of free 
energy, because it does not resolve the monomer con- 
stituents of droplets. Instead, one can fix the origin of 
free energy profiles in the 2d Ising model by requiring 
that Eq. Q returns the free energy of (say) clusters of 
size 1. The latter quantity can be calculated exactly in 
the Ising model, so fixing d(T) [TO] . 



In Fig. [T] we compare free energy profiles computed 
from umbrella sampling simulations with Eq. Pan- 
els (a) and (c) show droplet free energy profiles for var- 
ious values of k^T, for magnetic field strengths ft = 0.1 
and ft = 0.06. Open circles denote the result of um- 
brella sampling simulations; solid lines are obtained us- 
ing Eq. ([3| with the Shneidman et al. effective surface 
tension. As per Ref. [10 , the theoretical prediction fits 
the simulation data well across the range of conditions 
studied. Panels (b) and (d) compare free energy bar- 
riers obtained from simulation and theory in the plots 
directly above them. Here the dotted lines indicate the 
barrier height predictions of the uncorrected CNT ex- 
pression, Eq. ([3| with d = r = 0. In all cases the un- 
corrected theoretical prediction underestimates the bar- 
rier height. Upward pointing triangles show the bar- 
rier height prediction of Eq. Q using the Onsager so- 
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FIG. 2: Free energy profiles for growing a droplet in the 
presence of a wall, in the 2d Ising model (lattice gas repre- 
sentation), for various particle- wall attraction strengths J s . 
Nucleation in the absence of a particle-wall attraction is as 
in the bulk: curves for J s = 0, 0.2 and 0.4 coincide with 
the bulk curve (dashed red line) . The nucleation barrier and 
critical nucleus are reduced for large enough particle-wall 
attraction. The arrow marks the lattice gas particle- wall at- 
traction (J s = J/2) that is equivalent to the Ising model rep- 
resentation in the presence of an inert wall (K s = 0; see main 
text for discussion). Snapshots above are representative of 
critical nuclei at the particle- wall interactions specified. 



lution for the surface tension. Agreement between it 
and simulation is good at high temperature, but not at 
low temperature, where droplets are anisotropic. Also 
shown are averaged droplet profiles for critical nuclei at 
the extreme values of k^T at each value of h. These 
were obtained by averaging over all configurations for 
which N = N c during umbrella sampling. Snapshots 
are scaled relative to the size of the largest cluster shown 
(which occurs when k^,T = 0.9 and h = 0.06). Droplets 
are noticeably non-circular at low temperatures. 



B. Nucleation at a planar surface is faster than in 
the bulk only if the surface is sufficiently attractive 

With con fidence in our sampling protocol (Ap- 
pendix 



A 1 ) established, we next turn to the ques- 
tion of how a planar surface affects nucleation. Fig. [2] 
shows free energy barriers to nucleation for the 2d lat- 
tice gas (Eq. ([I])) in the presence of a flat wall, for a 
range of particle- wall interaction strengths J s . We set 
J = 3.2, ii = -6.3 (equivalent to K = 0.8, h = 0.05 
in the Ising model representation) and k^T = 1. Under 
these conditions the particle (up spin) phase is ther- 
modynamically preferred to the initial vacancy (down 



spin) phase. In the absence of a particle- wall attrac- 
tion, particles are effectively repelled by the wall, for 
reasons of entropy (sites available in the bulk exceed 
those available near the wall) and geometry (the nu- 
cleus shape that minimizes the surface-to-area ratio in 
two dimensions, a circle, can form only in the bulk). For 
particle-wall attractions not strong enough to overcome 
the entropic penalty of wall confinement, nuclei again 
grow in the bulk of the simulation box. (Sampling was 
initialized using wall-hugging droplets generated using 
large J s ; when bulk nucleation was preferred, droplets 
moved away from the wall). For sufficiently large at- 
tractions nuclei do grow at the wall, and the free energy 
barrier to nucleation and the size of the critical nucleus 
are smaller than their bulk counterparts. A substantial 
particle-wall attraction is needed to counter the favor- 
able entropy associated with bulk nucleation: in other 
words, a surface will enhance nucleation only if it pos- 
sesses a sufficiently strong attraction for the nucleating 
phase. 

While this observation is intuitively reasonable, we 
note that it is much more apparent in the lattice gas 
representation than the Ising one. The authors of 
Refs. [9j [32] studied nucleation using the Ising Hamil- 
tonian Eq. Q augmented by a bulk-wall interaction 

-\wall 



K s SiSJ, for the particular case K s = 



T^iwaH 
Ising 

(i.e. an energetically inert wall). Although K s = 
means that the wall has no energetic preference for ei- 
ther phase, the repulsion between unlike spins in the 
bulk leads to an effective attraction between the nucle- 
ating phase (in those papers the up-spin phase) and the 
wall [45 . Carrying through the Ising-lattice gas trans- 
formation, it can be shown that an Ising model in con- 
tact with an energetically inert wall (K s = 0) is equiva- 
lent to a lattice gas in contact with a wall that possesses 
a substantial interaction for particles, i.e. Eq. ([I]) with 
J s = J/2 [33]. This limit is marked by an arrow in 
Fig. [2] In both cases, one must engineer a substantial 
attraction between the nucleating phase and a wall be- 
fore nucleation happens at the wall in preference to in 
bulk (a result confirmed by simple scaling arguments: 
see Appendix [B|. In the lattice gas representation the 
coupling J s might be regarded as a literal wall-particle 
attraction; in Ising language, the attraction between up 
spins and the wall can be regarded as an effective one, 
mediated by 'solvent' (the down-spin phase). 



C. Nucleation in 2d pores can be faster still 

If a surface is sufficiently attractive, then, it can ren- 
der nucleation faster than in the bulk. Nucleation in 
a pore made out of that surface can be faster still, be- 
cause pore corners provide a convenient initiation site 
for the new thermodynamic phase [9]. Further, for given 
thermodynamic conditions, Ref. [9] demonstrated that 
there exists a pore size that maximizes nucleation rate. 
The existence of this maximum follows from the fact 
that as one makes a pore bigger, the rate for nucleation 
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FIG. 3: (a) Free energy profiles for nucleation in 2d Ising model pores. Profiles are consistent with the two-step mechanism 
described in Ref. [9], revealing a barrier to nucleation into the pore, and a barrier to nucleation from a filled pore into the 
bulk (see snapshots). The width w of the pore governs the heights of the two barriers, which show opposing dependencies 
on w. The total nucleation rate is a competition between these two processes and is optimized, for given thermodynamic 
conditions, by a specific pore width. The expanded region of the blue (w = 5) free energy curve illustrates that completed 
rows represent local metastable minima during the post-critical filling of a pore, (b) Barrier to nucleation inside a pore 
(black triangles) and out of filled pore (red circles) as a function of the width of the pore w. Note that the pore width that 
maximizes nucleation rate cannot be determined directly from the intersection of the two curves; instead, one must compute 
nucleation rates explicitly. The inset shows nucleation rates, i?i n ,out, computed using the forward flux sampling method [21]. 
In agreement with the results of Ref. 9 , we find that the overall nucleation rate (the reciprocal of the sum of the 'in' and 
'out' nucleation timescales) is optimized for a pore about 12 sites wide. 



into the pore is reduced, while the rate for nucleation 
out of the pore (into solution) is enhanced. In Fig. |3ja) 
we show free energy profiles for nucleation in a 2d pore. 
These profiles complement the nucleation rate calcula- 
tions of Ref. [9 j, confirming the existence of a double 
barrier to pore-mediated nucleation into solution. Un- 
der our umbrella sampling protocol, nucleation first oc- 
curs within the pore (starting in one of the corners due 
to the greater number of favorable energetic contacts 




x 



FIG. 4: Geometry for CNT-like scaling argument for pore 
nucleation 9 : (a) quarter-circular droplet nucleating in a 
pore corner, and (b) semicircular droplet nucleating out of a 
filled pore. Simple approximations for the free energy barri- 
ers in cases (a) and (b) show that the barrier for nucleation 
within a pore increases with w, while the barrier for nu- 
cleation out of a pore decreases with w. This competition 
implies the existence of an optimal pore with. 



there), and is followed by nucleation from a filled pore 
into the bulk. As highlighted in the boxed region, free 
energy profiles for narrow pores show local metastable 
minima during the post-critical filling of a pore. Each 
minimum corresponds to a filled row of the pore. 

As described in Ref. [9 j, the barrier to nucleation 
within the pore increases with increasing pore width, 
while the barrier to nucleation out into solution shows 
the opposite trend (see panel (b)). It should be noted 
that the pore width that maximizes nucleation rate can- 
not be determined directly from the intersection of the 
two curves: one must compute droplet nucleation rates 
explicitly. The inset shows the two nucleation rates, 
^in,out? m units of Monte Carlo steps per bulk lattice 
site, computed using forward flux sampling [21 j. In 
agreement with the results of Ref. [9 we find that the 
overall nucleation rate (the reciprocal of the sum of the 
'in' and 'out' pore nucleation timescales) is largest for a 
pore approximately 12 sites wide. 

A simple CNT-like approximation confirms that, 
given a sufficiently strong particle- wall attraction, there 
must exist a pore size that maximizes nucleation rate 
(this argument is that of Ref. [9 , modified to account 
for variable particle- wall surface tension). Approximat- 
ing the in-pore nucleus as a quarter-circle droplet of 
radius x growing from the corner of a pore of width w 
(see Fig. [4^a)) suggests a free energy cost for the nu- 
cleus of G- in (x) c± 2xcr s + 7rcrx/2 — Ag£ 2 7r/4 (here <r s 
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FIG. 5: Heterogeneous nucleation in three dimensions in the presence of a cuboidal pore of L z — 10 lattice sites and aspect 
ratio L x : L y . (a) Contour map showing the barrier to nucleation inside a pore. Within the white region there is no stable 
nucleus within the pore, (b) Contour map showing barrier to nucleation out of a pore. Within the white region there is no 
barrier to nucleation out of the pore, i.e. a critical nucleus can form within the pore and grow without bound, (c) Contour 
map of the largest barrier to nucleation (either in pore, out of pore, or single barrier). Since nucleation occurs in a pore 
corner, and the number of corners does not change with increasing pore size, the nucleation barrier is strongly dependent 
upon only a single horizontal dimension. If one pore dimension takes this 'correct' size, the nucleation barrier depends only 
weakly upon the other pore dimension, provided that the latter is large enough. The region of optimum pore geometry 
therefore forms a well-defined band, as highlighted by the dashed white line. The color scale is the same for all plots. Panels 
(d)-(f) show complementary contour maps of the nucleation rate (calculated using forward flux sampling) for the panel 
directly above. The low lying band seen in the panel (c) corresponds, to within a lattice site, to a ridge in panel (f), along 
which the overall rate of nucleation is maximized. 



is the droplet- wall surface tension). Assuming that the 
pore is narrow (so that the function G- m (x) does not 
reach its turning point for x < w) then the barrier to 
nucleation goes as G™*(w) ~ (2a s + ttg/2)w - G(w 2 ), 
which increases (sub-linearly) with pore width w. For 
nucleation out of a filled pore (Fig. Qb)), by a semi- 
circular droplet radius the free energy profile can be 
approximated as G out (x) ~ (2x— w)a s +7rxa — Agx 2 7r/2, 
which gives a pore width-dependent nucleation barrier 
of G™ a t x (w) = (2cr s + 7rcr) 2 /(7rA^)-w. This decreases 
linearly with w (note that the 'out' barriers seen in Fig. [3] 
are indeed approximately linear in w). An optimum 
pore width arises naturally from the competition be- 
tween these two processes. 

We stress, however, that the degree of attenuation of 
the nucleation barrier due to the pore depends on the 



particle-pore attraction, and can range from nothing at 
all (for small J s ), to total (for large J s ). 



IV. RESULTS: 3D SIMULATIONS 

A. Nucleation barriers in 3d pores depend on 
pore size and aspect ratio 

The analysis of Ref. [9] can be straightforwardly ex- 
tended to treat cuboidal pores. We studied heteroge- 
neous nucleation in a three-dimensional lattice gas of 
30 3 sites, for J = 1.6, /i = —3.5, in the presence of a 
cuboidal pore of dimensions L x x L y x L z embedded 
in a surface. Both the pore and the surface possessed 
stickiness parameter J s = 0.6. For a fixed pore depth 
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FIG. 6: (a) Largest barrier height and total nucleation rate taken from diagonal cuts (L x — L y ) along contour maps (c) and 
(f ) in Fig. [5] The square aspect ratio pore that maximizes the nucleation rate has a side of length L x — L y — 6. (b) If L x 
takes a certain, optimal value, then the barrier height depends only weakly on L y , provided that L y is large enough, (c) The 
same is true of the overall nucleation rate. 



of L z = 10 sites we computed nucleation free energy 
profiles. Fig. [5] demonstrates that, in general, there ex- 
ists a barrier to filling a pore (Gin, Fig. [5|a)), and for 
growing a droplet from a filled pore into solution (G ut> 
Fig.gb)). For very large pores, however, droplets can 
attain criticality in the pore, and no barrier to growing 
into solution exists. Very small pores, by contrast, can- 
not be filled without substantial free energy cost, and 
offer no enhancement of nucleation relative to a planar 
surface. When such pores are present, nuclei grow on 
the surrounding flat substrate (made of the same mate- 
rial as the pore). 

Fig.[5jc) shows a contour plot of the larger of the 'in' 
and 'out' barriers, G max , as a function of a pore's hori- 
zontal dimensions. All pores provide a significant atten- 
uation of the free energy barrier to nucleation relative to 
bulk, because the pore and its surrounding surface is at- 
tractive (here the bulk nucleation barrier is ~ 67&bT). 
However, G max varies dramatically with pore geometry. 
Achieving maximum attenuation of the nucleation bar- 
rier requires making only one of the pore's horizontal 
dimensions the 'correct' size, provided that the second 
is large enough. Consequently, the region of optimal 
pore geometry forms a well-defined band on the con- 
tour plot. Panels (d)-(f) show complementary contour 
plots of nucleation rates, computed using forward flux 
sampling. The band seen in panel (c) corresponds, to 
within a lattice site, to a ridge in panel (f), along which 
the overall rate of nucleation is maximized. 

A selection of cuts along the contour plots (corre- 
sponding to varying either the size or shape of a pore) 
are shown in Fig. [6j As is evident from panels (b) and 
(c), the nucleation barrier and rate are not strongly de- 
pendent on the larger of a pore's dimensions, as long 
as that dimension is large enough. This can be under- 
stood by noting that nucleation within a pore occurs 
at a corner; as a droplet grows, it is stabilized energeti- 



cally when it encounters the closer of the two other pore 
walls. Provided that the distance to the further wall is 
large enough, that distance does not strongly affect the 
ability of the nucleus to grow out into solution. 

This observation suggests that it is preferable to 
pattern a surface by repeating in it copies of a well- 
chosen pore (a small pore that effects a substantial re- 
duction in nucleation free energy barrier, such as the 
10 x 3 pore), rather than e.g. etching long grooves in 
it. Consider, as one possible choice, the best square 
pore (L x = L y = 5; see Fig. [6|a)). The larger of 
the in- and out barriers to nucleation for that pore is 
« 20/cbT. For a substrate row of (large) length L, 
we therefore expect the nucleation timescale associated 
with repeated, closely-spaced copies of the square to 
be T square ~ (6/L) exp(20). By comparison, the barrier 
to nucleation for a periodic groove of the same width, 
built by imposing periodic boundaries in the ^-direction 
of a box of length 30, is 24.7 k^T, and so we expect 
the nucleation timescale for a groove of length L to be 

'groove rsJ 

(30/L)exp(24.7). We therefore expect that 
replacing a single long groove of width 5 by an array of 
square pores of the same width will increase nucleation 
rate by a factor of r groove /r sqiiare ~ 500. Rate calcu- 
lations done using forward flux sampling are consistent 
with this estimate. 



B. Nucleation mechanisms in 3d pores also 
depend on pore size and aspect ratio 

Fig. [7] illustrates the range of behaviors associated 
with different pore shapes and sizes. For some pores, 
nucleation happens in a single step: 1) droplets are not 
stable within very small pores; instead, a critical nu- 
cleus appears on the surrounding surface (top left); or 
2) droplets can attain criticality within very large pores 



8 




FIG. 7: Nucleation mechanism as a function of pore aspect ratio (L x ,L y ) for pores of fixed depth L z = 10. Snapshots 
and free energy profiles illustrate the key behaviors seen: (top left) droplets are unstable within the pore, and so a critical 
nucleus can only form on the flat surface instead; (top right) the pore is so large that a droplet within it becomes critical 
before the pore is filled; (bottom left) the filled pore is metastable with respect to the empty lattice (solution) and full lattice 
(stable, nucleated phase); (bottom right) the filled pore is stable with respect to solution and metastable with respect to the 
nucleated phase. For clarity, snapshots show only particles within the largest cluster in the simulation box. 



(top right). For other pores, nucleation happens in two 
steps: 3) small filled pores are metastable with respect 
to both solution and the nucleated phase (bottom left); 
and 4) filled, moderately-sized pores are stable with re- 
spect to solution, and metastable with respect to the 
nucleated phase (bottom right). 



V. CONCLUSIONS 

We have studied homogeneous and heterogeneous nu- 
cleation in the 2d and 3d Ising models. Our key result 
is the extension of the work of Ref. [9 to calculate rates 
and free energy profiles for nucleation in the 3d Ising 
model in the presence of cuboidal pores. Pores of well- 
chosen aspect ratio can dramatically speed nucleation 
relative to a planar surface made of the same material, 
while badly-chosen pores provide no such enhancement. 
Further, for given thermodynamic conditions, and a suf- 
ficiently strong pore-particle attraction, there exists a 



pore size and aspect ratio ideal for promoting nucle- 
ation. Fig. [5] summarizes the importance of pore choice 
in reducing free energy barriers to nucleation. A suf- 
ficiently attractive surface can dramatically reduce the 
nucleation barrier relative to that in bulk. For the pa- 
rameters considered here, the bulk free energy nucle- 
ation barrier is about 67/cbT, while the barrier in the 
presence of an attractive wall ( J s = 0.6) is about 45 k^T. 
A pore made from the same material can reduce the 
barrier even further: the barrier is about 20 k^T for 
the optimally-sized square pore (Fig. j8^a)). However, a 
badly-chosen pore offers no improvement over a planar 
surface: although small droplets appear first within the 
pore of Fig. [8^b), the critical nucleus forms instead on 
the surface surrounding it. 
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FIG. 8: Well-chosen pores can dramatically quicken nucleation relative to a planar surface made of the same material, while 
badly-chosen pores do little except take up space on the substrate, (a) If a planar surface has a sufficient attraction for the 
nucleating phase then it can promote nucleation relative to the bulk case. Here the bulk nucleation barrier is about 67A;bT; 
at an attractive wall (J s = 0.6) the barrier is instead 45/cbT. A pore of well-chosen geometry can reduce the barrier even 
further. Here nucleation within a square 5x5 pore has a maximum barrier of about 20 k&T. (b) By contrast, a badly-chosen 
pore offers no enhancement of nucleation over a flat surface of the same material. Given a pore so narrow that droplets are 
unstable within it, nucleation takes place instead on the surrounding surface, (c) Summary: nucleation at a planar surface 
(of stickiness J s — 0.6) is 9 orders of magnitude faster than in the bulk. A long groove of optimum width promotes nucleation 
rate by a further 7 orders of magnitude. A square pore of optimum width gives rise to nucleation that is about 70 times 
faster still, suggesting that, for given thermodynamic conditions, a raster arrangement of pores - repeating copies of the 
square - is a better way to speed nucleation than e.g. scoring long grooves in the surface. 
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Appendix A: Umbrella sampling 

1. Sampling 

Several implementations of umbrella sampling [20] 
for the study of nucleation are described in the liter- 
ature [35ti37| . We used a hybrid of the following two 
methods. In the first method (Method A) we mea- 
sured the distribution of sizes of all connected clusters 
in the simulation box [T6l [34j [35] - We carried out 'hard 
wall' umbrella sampling simulations [35 , constraining 
the simulation to a 'window' (of length 10 and with an 
overlap of 5 with its next neighbor) by rejecting any 
spin flip that made the largest cluster in the simulation 
box larger or smaller than the window's limits. Within 
each window we recorded histograms of the density of 



clusters of all sizes that fall within the window's bounds, 
measuring p(N) = (Mn)/V, where Mjy is the number 
of clusters of size TV in the simulation box, and V is 
its volume. We combined data from all windows using 
the weighted histogram analysis method (WHAM) [38 , 
giving the free energy G&(N) = — kBTlnp(N). We 
checked this sampling procedure by reproducing, for the 
3d Ising model, the curves shown in Fig. 5 of Ref. [35] 
and Fig. 1 of Ref. [16]. 

While simple to implement, the hard wall umbrella 
sampling technique can become inefficient when there 
exist steep gradients in free energy. In this situation 
the sampling within a window is largely confined to the 
region adjacent to the wall lower in free energy, leav- 
ing the other end of the window badly sampled. Al- 
though it is possible to circumvent this problem by in- 
creasing the number of sampling windows (or their de- 
gree of overlap), this increases the computational cost 
of the method. A common solution (and the second 
method considered, Method B) is to constrain the size 
N max of the largest cluster using a harmonic bias po- 
tential, ki(N max - N{ arget ) 2 /2 [36]. Here i = 1,2,... 
designate the different sampling windows, ki is a spring 
constant, and A^ arget is the cluster size at the center 
of window i. The width of a window is determined 
by the strength of the spring constant ki which can be 
softened or stiffened according to the local free energy 
gradient in order to improve sampling. We generated 
histograms P B (^max) = A/"(7V max )/ E Nmax ^(7V max ) by 
adding 1 to a register J\f(N max ) if, after every trial move, 
the largest cluster in the system was of size iV max (re- 
gardless of how many clusters of that size there were). 
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FIG. 9: Artifacts of measuring only the size of the largest cluster in the system are pronounced at deep supercooling but not 
at shallow supercooling. In panel (a) we show free energy profiles calculated by umbrella sampling methods A (measuring 
sizes of all clusters) and B (using as a reaction coordinate the system's largest cluster), under conditions that give rise to the 
smallest barrier in Fig[TJa). Here artifacts are apparent in the small- N data produced by method B, because it ignores the 
many small clusters present in the simulation box (the inset, generated using unconstrained simulations, shows the likelihood 
that the largest cluster in the system is of size N). For larger N the shape of the curves generated by the two methods agree, 
and when the 'B' curve is shifted up by about InV [34], it sits on top of the 'A' curve, (b) For conditions that give rise to 
the largest barrier in Fig.JTJb), by contrast, the curves generated by the two methods are similar, because the average cluster 
size is 1. 



To generate overlapping histograms we used a window 
spacing of ^target — iV t * arget = 5 and a spring constant of 
ki = 0.2. Sampling was done for a minimum of 10 6 MC 
sweeps within each window. An MC sweep consisted of 
^Vbuik attempted trial moves, where iVbuik was the num- 
ber of lattice sites in the bulk of the simulation box, i.e. 
sites that are not part of a wall or pore. We updated 
and recorded cluster size distributions after every trial 
move, so maximizing the efficiency of our simulations. 
We used umbrella integration [39 to unbias the results 
of umbrella sampling simulations, and used WHAM [38] 
to check this procedure and resolv e fine details of certain 
free energy curves (see Appendix | A 2[ ). This procedure 
gives the 'free energy' G B (N max ) = -feTlnP B (^Vmax). 

Unlike method A, which measures the distribution 
of sizes for all clusters, method B instead samples the 
probability that the largest cluster is of size iV max . As 
pointed out by Maibaum [34., a free energy penalty is 
incurred whenever Af max is constrained to sizes smaller 
than the average cluster size, iV av , seen in unconstrained 
simulations (that do not result in nucleation). Fig. [9] 
shows a comparison between the two umbrella sampling 
methods at conditions of deep (a) and shallow (b) su- 
percooling. Here AG is defined as G\(N) — Ga(1) and 
GB(^V max ) — min(GB (N max )) for methods A and B re- 
spectively. At deep supercooling the average cluster size 
at small N is 5; consequently, a spurious increase in AG 
is seen in the free energy curve obtained using method 
B for N < 5 [34 . In contrast, at shallow supercooling 
the average cluster size at small N is 1, and the methods 
agree. 

Although small-TV artifacts can be present using 



method B, it is important to note that the shape of 
the free energy curve is correct for N ^> N av . We can 
therefore take advantage of the improved sampling that 
method B provides by using it to sample large clus- 
ters, and using method A to gather statistics for small 
clusters (which is typically cheap to do). The curves 
reported in the text were stitched together using both 
methods, with the method B results shifted vertically 
to match the small-cluster data obtained using method 
A. Comparison with the CNT-like predictions in Sec- 
tion III A and the results of Ref. [10] allow us to confirm 
that this scheme works. 

To complement our free energy sampling we calcu- 
lated nucleation rates directly using the forward flux 
sampling (FFS) method [2T] , FFS is reasonably insen- 
sitive to the choice of reaction coordinate [21 , and we 
verified that consistent results were obtained using as a 
reaction coordinate 1) the size of the largest cluster in 
our simulation box and 2) the total number of particles. 
All interfaces were spaced 10 particles apart, and 10000 
crossings were stored at each. Rates were measured in 
units of Monte Carlo steps per bulk site. 

To measure nucleation rates in the presence of pores 
we followed the protocol outlined in Ref. [9] and decom- 
posed the overall nucleation rate into two parts: the rate 
for pore filling i? in ; and the rate for nucleation out of 
an already filled pore R out . These rates define the mean 
timescales for the two nucleation processes, r- in = iir 1 
and r out = ^out? which can be combined to give the 
overall nucleation time, r to tai = + Tout, and hence 
the overall nucleation rate i?totai = r totai- 
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2. Unbiasing distributions 



Appendix B: Nucleation barriers at a planar 
surface are reasonably well described by CNT 



We combined data from each window using the 
weighted histogram analysis method (WHAM) [38] 
and/or umbrella integration [39]. WHAM takes as its 
input an overlapping sequence of probability distribu- 
tions, and effects a self-consistent iteration to reconsti- 
tute the underlying free energy curve. Umbrella integra- 
tion, by contrast, assumes Gaussian probability distri- 
butions within each window, and takes as input only the 
mean and variance of the sampled distribution. It also 
does not require windows to overlap. It is therefore com- 
putationally cheaper to implement than WHAM. For 
example, the contour plots shown in Fig. [5] were made 
from 210 sets of simulations, each comprising 200 indi- 
vidual sampling windows. Using umbrella integration 
it was possible to compute all of the 210 free energy 
curves in less than 10 seconds. By contrast, WHAM 
took north of 20 minutes to process an individual free 
energy curve. 

As the authors of the method caution [39 , umbrella 
integration misses features of a free energy curve if the 
sampling within windows is non-Gaussian, as happens 
when a pore fills layer by layer. In such cases, WHAM 
reveals subtle local features (metastable minima) in the 
free energy profile: see e.g. the w — 5 profile in Fig.[3ja), 
enlarged in Fig. [To] for clarity. Umbrella integration 
does not. However, umbrella integration gets correct 
the overall shape of the free energy profile, including 
the positions and heights of the large free energy barriers 
(corresponding to the initial nucleation event within the 
pore, and nucleation of a droplet out of the filled pore). 
We therefore used umbrella integration to broadly sur- 
vey free energy landscapes in Figs[5}{7| and used WHAM 
to resolve subtle features of particular curves. 



The requirement that a wall attract a nucleating 
phase before bulk nucleation is suppressed is consistent 
with a classical nucleation theory-like scaling argument. 
We can estimate the barrier for nucleation at a wall by 
considering a droplet whose shape is the portion of a 
circle of radius R that lies above a surface when the 
surface-circle intersections subtend an angle i/j at the 
center of the circle (see Fig. 11 1 inset). When ip = the 



circle just touches the wall and we have nucleation in 
the bulk; when ip = tt the droplet is a semi-circle; when 
ip = 2tt we consider the droplet to wet the wall. We 
estimate the free energy of the droplet as 



G s (R,ip) = -AAg + al- 



cT s l s ^-e^)\nN ho ^ (Bl) 



The first three terms account for droplet area A = N = 
(tt — ip/2)R 2 + (1/2)R 2 sin^; curved droplet perimeter 
I = (27r — iji)R\ and wall-contacting droplet perimeter 
l s = 2Rsm(ip/2). <j s is the wall-droplet surface tension. 
The final term = 1 if ij) > 0, and is zero if ijj = 0) 

accounts for the fact that there are more ways of placing 
a droplet in the bulk than at the wall of the system, 
i.e. there is an entropic cost associated with movi ng a 
cluster from the bulk to the wall. Maximizing Eq. (Bl) 
with respect to R gives the critical radius 



a (2tt - + 2a s sin(^/2) 
A#(27r-^ + sin(V0) 



(B2) 



Substituting this expression into Eq. (Bl) gives the free 
energy barrier to nucleation: 
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(B3) 



Upon setting i/j = 0, we recover the conventional bulk 
solution, Eq. ([3| with r = d = 0. The droplet makes a 
contact angle with the surface of = tt — Balanc- 
ing surface tensions using Young's equati on [40] gives 
a cosO 



-<r s . Using these results in Eq. (B3) gives a 
contact angle-dependent free energy barrier 



Gi 



where 



bulk 



(i? bu lk)+ l e(7r . 



0)lniV box , (B4) 



FIG. 10: An enlargement of the boxed region shown 
in Fig. |3ja). We used umbrella integration as a fast way to 
reconstitute the large-scale features of free energy profiles, 
and WHAM to resolve profiles' small features when proba- 
bility distributions within windows were non-Gaussian. 



f{0) = -[e- -sin(20) . 



(B5) 



This is the Turnbull estimate for the free energy of a 2d 
droplet at a surface [41 j. In terms of surface tensions 
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the function / reads 



f(a,a a ) 



l r<r ; 



;y/{(T- cr s )(cr + a s ) + arccos(-a s /cr) . 

(B6) 

To compare this estimate with the results of simula- 



tion we added to Eq. (B4) the difference between the 



numerical value for the bulk free energy barrier, calcu- 
lated via umbrella sampling, and the uncorrected CNT 
prediction, i.e. 



Gi 



-+ f(0)G hulk (R h c ulk ) + -e(7r-0)\nN ho 



-AG 



bulk, 

(B7) 

where AG bulk = GftSk - G^ lk r | d=r=0 . This modifica- 
tion assumes that the ill-defined origin of our CNT-like 
expression can be fixed by requiring that in the bulk 
limit it return s the results of computer simulations (see 
Section III A). This modification is ad- hoc and uncon- 
trolled. 

We take Onsager's solution for the Ising model sur- 
face tension, Eq. Q, as an approximation [46] for the 
droplet-solution surface tension: 



a = J/2- k^T In [coth(/3J/4)] . 



(B8) 



Note that J is the lattice gas coupling, not the Ising 
one. To estimate a value for the droplet-wall surface 
tension, we note that the two terms in the Onsager sur- 
face tension account for the energy per unit length of a 
planar interface (first term) , and the free energy per unit 
length of fluctuations normal to that interface (second 
term) [42]. To see this, consider an interface of horizon- 
tal length L between up- and down spins in a 2d Ising 
model (at h = 0), where the vertical position of the in- 
terface at the n th lattice site across it (n = 1, 2, . . . , L) 
is u n . The energy cost of fluctuations normal to the 
interface is therefore 2ifJ^ =1 |iz n +i — u n \, and the as- 
sociated partition function is 



^interface — ^ ^ ^ 
{u n } 

/ oo 



2^E n \u n + l-U n \ 



-2K/3\Au\ 



\u= — oo 



coth L {(3K) . 



(B9) 



The free energy per unit length associated with per- 
pendicular fluctuations of the interface, /interface = 



-L-^eTln 

^interface, is therefore the second term in 
Eq. ([B8| (recall that K = J/4). We now guess that 
when the wall is attractive enough that a droplet re- 
mains in close contact with it, the wall-droplet interface 
acts as if the wall can fluctuate. Clearly this is not true 
microscopically, and is likely to be a poor approximation 
if the droplet surface moves appreciable away from the 
wall. We nonetheless conjecture that the surface tension 
between the droplet and the wall can be approximated 
as 

a s = ( J - J s )/2 - /c B Tlncoth(/3(J - J s )/4). (BIO) 



We shall use Eqns. (B8 



and ( |B1Q ) to relate the surface 



tensions entering Eq. (B6) to the particle- wall attrac- 



tion J s used in our simulations. Finally, our CNT-like 
prediction for the barrier to nucleation in the presence 
of a surface is 



G max — min (G^, G^ k ) 



(Bll) 



i.e. we take the bulk barrier whenever the barrier to 
nucleation at the surface is larger than it (by virtue of 
the entropic penalty of wall confinement). Fig. [TTj^a) 
shows a comparison between theoretical (solid linejand 
simulated (circles) free energy barriers. Both calcula- 
tions show that nucleation only occurs at the wall when 
the particle-wall attraction is strong enough to offset 
the entropic penalty of removing the droplet from the 
bulk. Moreover, the quantitative agreement between 
the two methods is reasonable, which is surprising in 
light of the crude nature of the approximations we have 
made. Fig. [TTj^b) shows the umbrella sampling free en- 
ergy barrier as a function of the particle- wall interaction 
strength. Fig. [Tljc) shows the contact angle computed 
from Young's equation, which is 180° whenever bulk 
nucleation is preferred. The numbered snapshots show 
the average profiles of critical nuclei at different values 
of the reduced surface tension. Those on the left are 
generated using the contact angle from the solution to 
Young's equation, while those on the right are averages 
from umbrella sampling simulations. In all cases the 
agreement between theory and simulation is good. Here 
we have made a pictorial comparison between droplet 
shapes predicted by theory and computed by simula- 
tion; previous work demonstrates quantitatively that 
CNT can predict the Ising model droplet contact an- 
gle 03] El]. 
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